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We develop a non-local density functional by a direct modeling of the shape of exchange-correlation 
(xc) hole in inhomogeneous systems. The functional is aimed at giving an accurate xc-energy and an 
accurate corresponding xc-potential even in difficult near-degeneracy situations such as molecular 
bond breaking. In particular we demand that: (1) the xc hole properly contains —1 electron, (2) 
the xc-potential has the asymptotic — 1/r behavior outside finite systems and (3) the xc-potential 
has the correct step structure related to the derivative discontinuities of the xc-energy functional. 
These demands are achieved by screening the exchange hole in such a way that the pair-correlation 
function is symmetric and satisfies the sum-rule. These two features immediately imply (1) and (2) 
while the explicit dependence of the exchange hole on the Kohn-Sham orbitals implies (3). For the 
dissociating hydrogen molecule, we obtain much improved xc holes and thus also improved binding 
energies at different interatomic distances, as compared to the local density approximation. 



I. INTRODUCTION 

The local density approximation (LDA) is the simplest 
functional in density functional theory (DFT) and has 
been around since the advent of DFT [1, 2]. Although the 
LDA has only a local dependence on the density, it has 
been tremendously successful in describing the ground- 
state properties of solids, surfaces and large molecules. 
Its shortcomings have been obvious from the beginning 
and an enormous effort has gone into the search for bet- 
ter approximations. This task has proved to be exceed- 
ingly difficult. One should, however, not be surprised 
or disappointed. Density functional theory provides a 
way to reduce the full, interacting many-body problem 
to a simple non-interacting one. Therefore, an accurate 
density functional for the total energy would provide a 
surprisingly simple way to solve the complicated many- 
body problem — at least as far as static properties are 
concerned. Nevertheless, in the past decades consider- 
able progress has been made. With the advent of the 
generalized gradient approximations (GGAs) [3, 4], bond 
lengths and atomization energies were greatly improved 
as compared to those of the LDA. Unfortunately, the 
exchange-correlation (xc) potentials of the GGAs have 
several undesirable features. In particular, they decay 
exponentially outside finite systems like the xc potential 
of the LDA and unlike a correct —1/r decay. Conse- 
quently, neither the LDA nor the GGAs produce proper 
Rydberg levels [5] and ionization potentials are too low. 

An important class of extensions to the GGAs came 
from the observation that exchange energies are much 
larger than correlation energies. This suggests that ex- 
change should be treated exactly, while using a GGA 
only for the correlation energy. Full inclusion of "exact 
exchange" does however not work well in practice, since 



there is a large cancellation of errors between the LDA 
or the GGA versions of the exchange and correlation en- 
ergies. On the other hand, using only a portion of exact 
exchange combined with a GGA does yield quite accurate 
bond energies in molecules [6]. An important advantage 
of "exact exchange" is that it provides some necessary 
improvements of the xc potential. For instance, the in- 
clusion of exact exchange gives a stepped structure in 
the xc potential [7, 8] thus improving the description of 
the atomic shells as well as enabling a neutral dissocia- 
tion of heteronnclear molecules [9-11]. It also gives an 
xc-potential with the proper asymptotic (—1/r) behavior 
which gives rise to Rydberg levels and a highest occupied 
KS eigenvalue in better agreement with the negative of 
the ionization potential. Usually, however, only a frac- 
tion of full exchange is incorporated in the so called hy- 
brid functionals meaning that the desirable features men- 
tioned above are only partially obtained. It seems that 
only functionals with a massive amount of fitting param- 
eters are able to handle 100% "exact exchange" [12, 13]. 
We mention here that inclusion of "exact exchange" is 
not mandatory in order to have good properties of the 
xc potential. The proper step structure as well as the 
correct asymptotics away from finite systems can also 
be obtained by modeling the potential directly [5, 14]. 
Such model potentials can indeed provide good response 
properties [15, 16] but they cannot easily be written as 
the functional derivative of some accurate functional for 
the xc energy. As a result, their implementations have 
been limited. Further improvements to the energy func- 
tional are also sought by adding the kinetic-energy den- 
sity to the functional arguments [17, 18] and including 
parts from many-body perturbation theory [19-23]. 

In this paper we would like to consider an alterna- 
tive approach by looking directly at electron correlations 
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in real space. In a correlated system we can consider 
the conditional probability of finding an electron at some 
point in space, when the position of another reference 
electron is given. The difference between this function 
and the unconditional probability (which is simply the 
density) is defined to be the xc hole and the knowledge 
of this function [24] is sufficient to calculate the xc-energy. 



The effect of exchange and correlation is to dig a hole 
in the density around each electron and this hole misses 
one electron. We can say that the task of a good xc 
functional is to provide an accurate description of the xc 
hole. The LDA and GGA assume that this hole has a 
spherical shape and an extent given by the Wigner-Seitz 
radius {Airr^n — 3, n being the electron density) at the 
reference electron. Unfortunately for the LDA, xc holes 
are definitely not spherical [25, 26]. An important class 
of functionals which aim to lift these restrictions is the so- 
called weighted density approximations (WDA) [27-30]. 
These functionals avoid the spherical xc hole by digging 
a hole out of the real density rather than in the density 
at the position of the reference electron. A nice feature 
of the WDA is that the asymptotic behavior of its xc 
potential has a Coulombic asymptotic decay instead of 
an exponential behavior as in the LDA and the GGA. 
An important symmetry of the pair-density (r(r',r) = 
r(r,r')) is, however, broken. This causes an additional 
factor 1/2 in the asymptotic decay of the xc potential, so 
that it decays too fast (as — f/(2r)) [31]. 



Here, we will advocate an approach in a similar spirit 
as the weighted density approximation [32]. However, 
we will take care not to destroy the symmetry of the 
pair-density and therefore, the xc potential will automat- 
ically have the correct asymptotic — l/r behavior. Fur- 
thermore, important information on the physics of the xc 
hole is provided by the dissociation of molecules. In par- 
ticular, a proper localization of the xc hole around the 
reference electron for a dissociated molecule is a chal- 
lenging task. The failure of current approximations to 
achieve this is reflected in their consistent inability to 
properly describe the breaking of chemical bonds. The 
most important aim of our new functional will therefore 
be a bold one: the new functional has to be able to de- 
scribe molecular dissociation. 



The paper is outlined as follows. First (Sec. II) we 
will give a more detailed discussion on the background to 
motivate the construction of the screened exchange (SX) 
functional. The actual construction of the SX functional 
is discussed in section HI. In section IV we show the 
results and finally conclude in section V. 



II. MOTIVATION 

A. Symmetry and asymptotics 

We will start from an exact expression [33-35] for the 
exchange-correlation energy 

where n(r) is the density. This expression gives the exact 
xc contribution to the interaction energy of the system, 
provided g is the exact pair-correlation function g of the 
system, defined as 

r(r,rO 
(7(r, r ) = — ; — r- 
■^^ ' ' n{v)n{v') 

Here the diagonal of the two-body reduced density matrix 
is defined as 

r(r,r')^5]r(m,rV) 

<T.(T' 

where il)'^{va) and '0(rcr) are the usual field operators. 
The xc energy, however, also contains a correlation con- 
tribution to the kinetic energy which is most conveniently 
included by integrating the interaction energy with re- 
spect to the strength A of the Coulomb interaction — 
while keeping the density fixed at the fully interaction 
one (A = 1) [32-35], i.e. 

5(r,r')= /dA5A(r,r'). 

The physical picture of representing the xc energy in this 
manner is that an electron does not interact with the full 
density, but depending on its position, r, it sees an effec- 
tive density n(r')5(r, r') with iV — 1 electrons. Therefore, 
the part 

p,e(r'|r)=n(r')(g(r,r')-l) (2) 

in the xc energy (1) exactly describes the density of minus 
one electron, a hole, which is reflected in the sum-rule of 
the xc hole 

y'dr'p.,(r'lr) = Jdr' n{r'){g{r,r') - l) = -1. (3) 

The local density approximation (LDA) proceeds by us- 
ing the xc hole from the homogeneous electron gas evalu- 
ated for the density at the position of the reference elec- 
tron, so the pair-correlation function is approximated as 
^(r, r') « g/i(lr — r'l; n{r)). Furthermore, if the distance 
]r — r'] is large the pair-correlation function g only dif- 
fers slightly from one, so it is quite reasonable to replace 
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n(r') by n(r) in the xc-energy (1). Combining both ap- 
proximations, we obtain the expression for the xc energy 
of the LDA 



E: 



LDA 



dr / dr ri(r) 



, n(r){5,(|r-r'|;n(r))-l} 



2 / dr/dr n{v) 



dr n(r)exc {n{v)), 



(4) 



where the function £xc('t-) is just the exchange-correlation 
energy per electron of the homogeneous electron gas. 
This expression for the xc energy of the LDA shows 
explicitly that its hole is approximated by a spheri- 
cal one. As mentioned above xc holes are not spheri- 
cal [25, 26]. The original weighted density approxima- 
tion (WDA) [27, 28] improves on the shape of the xc 
hole by not replacing n(r') by ?T,(r) in the xc-energy (1). 
Since the sum-rule (3) for the xc hole is no longer triv- 
ially satisfied, an effective density, n(r), is used as input 
into the pair-correlation function instead of the density 
at the reference position. The xc energy in the WDA is 
therefore 



w- . 1 /dr /dr' n(r)-(^^Hg^(l^-r1; n(r))-l}^ 



(5) 



where the effective density n(r) should be found by sat- 
isfying the sum-rule for the xc hole (3) at every point r 



dr'n(r'){.g,.(|r-r'|;n(r)) -1} 



-1. 



Unfortunately, it was found that the pair-correlation 
functional of the homogeneous electron gas is not a good 
approximation to the pair-correlation function of inho- 
mogeneous systems as molecules and surfaces. Using a 
more localized function for 5 — 1, the results were signif- 
icantly improved [29, 30, 36]. 

From the expression of the WDA xc energy (5) we 
immediately notice that the integral kernel is still asym- 
metric in the coordinates r and r' and therefore breaks 
the important symmetry g(r, r') = g(r',r). This is not 
so important for the value of the xc energy. But for the 
corresponding xc potential, we obtain 

WDA,„^ _ 1 /■ w n{r'){gh{\v-r'\;n{v))-l) 



^r"(r)=2 /dr' 



1 f^^, n(r'){g.(|r-r'|;n(r'))-l} 



|r - r' 

\ |dr'|dr"n(r')^^'"^ 



g^(|r'-r"|;n) 5n{v") 
5n{v") (5n(r) 



The first term actually decays as — l/(2r) as is obvious 
from the sum-rule. The long-range behavior of the other 
two terms is not so obvious, but in practice they turn 
out to decay exponentially [29, 30]. Therefore, the xc po- 
tential decays too slowly (as — l/(2r)) compared to the 
proper — 1/r decay [28]. The incorrect long-range be- 
havior of the potential is expected to have a significant 
effect on properties which highly depend on a proper xc 
potential such as the ionization energy and Rydberg ex- 
citations [7, 16, 37]. This failure can be attributed to the 
broken symmetry of the pair-correlation function. Con- 
sequently, one of the requirements of our new functional 
will be that it should satisfy the proper symmetry of the 
pair-correlation function, i.e. ^(r, r') = ^(r',r). 



B. Step structure from exchange 

As mentioned in the introduction, the steps in the KS 
potential are important for a proper description of the 
atomic shell structure [7, 8] and also the neutral disso- 
ciation of hetero-nuclear molecules [9, 10]. It has been 
shown that the necessary steps for the atomic shell struc- 
ture are already featured by the exchange energy [7]. 
However, the necessary step in the dissociation of hetero 
diatomic molecules is a less clear case, since the spin- 
symmetry has to be broken to provide the necessary lo- 
calization [38, 39]. 

Closely related, the exchange energy also has the nec- 
essary features for the integer discontinuity [7, 10], since 
the exchange term often changes radically when crossing 
an integer number of electrons due to the usual idempo- 
tency of the KS density matrix. The corresponding hole, 
the exchange hole, can be expressed in spin-integrated 
quantities as 



Px(rlrrof) 



1 l7^(r, rrcf) 

2 n(ri.of) 



(6) 



where the spin-integrated KS density matrix is defined 

as 

7s (r, I"') = ^7(rcr,r'cr) 

^^(vI/,l^t(rV)^(ra)lvI/,), 

cr 

with ^1 s as the KS wavefunction. The KS density matrix 
can alternatively be expressed directly in terms of the KS 
orbitals 0fc(r) as 
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:(r,r')=^nfe0,(r)0^(r'). 



(7) 



where rik are the occupation numbers, being simply or 
2 in the non-degenerate case. The exchange hole has the 
convenient property that it already satisfies the hole xc 
hole sum-rule (3). 

Since exchange already satisfies a number of important 
properties, it is often used as a starting point to model 



FIG. 1. The different holes of the £[2 molecule at Rh-h = 1-4 Bohr and the reference electron at 0.3 Bohr to the left of the 
right nucleus along the bond axis (rrcf ~ (0,0,0.4) Bohr). The positions of the nuclei are indicated by the blue lines and the 
position of the reference electron is indicated by the red line. The left panel shows the exchange hole, px(r|rrGf) = — |cr3(r)p, 
the middle panel shows the correlation hole, pc(r|rrcf), which provides a small correction to have the more localized real hole, 
pxc(r|rrcf). 




FIG. 2. Similar to the previous plots, but now for Rh-h = 5.0 Bohr. The reference electron is still at 0.3 Bohr to the left of the 
right nucleus along the bond axis (rrcf = (0,0,2.2) Bohr now). The exchange hole, px(r|rrcf), remains completely delocalized, 
so requires an equally large correction from the correlation hole, pc(rjrrcf), to obtain the real hole, pxc(rirrcf), localized around 
the reference electron. 



the full exchange-correlation effects. Traditionally, one 
adds a correction term, correlation, defined as the differ- 
ence between the exact exchange-correlation quantities 
and the ones with one exchange taken into account. For 
example the correlation hole is simply defined as 

/Oc(r|r,.of) EE Pxc(r|rrof) - Px(r|rrof)- 

Although this approach has some appeal, it turns out 
to be inconvenient in practice. Especially considering 
bond breaking situations as shown explicitly in the next 
section. 



C. Bond breaking 

In Fig. 1 we have plotted the different holes, px(r|ri.of), 
Pc(r|ri.of) and pxc(r|ri.cf) for the H2 molecule at equilib- 
rium distance i?H-H = Rc = 1-4 Bohr. Note that the 
quantities with correlation are at full coupling strength 
(A = 1) and not the integrated ones. Ideally we would 
like to have shown the integrated ones, but obtaining 
them is a rather difficult task. Since the effect of the ki- 
netic energy is not very large on the total energy [40] , we 
expect that the holes do not differ too much as well. The 



reference electron is fixed at 0.3 Bohr to the left from 
the right nucleus. The calculations holes were calculated 
from full configurations interaction (CI) results using the 
Is, 2s, 3s, 2p and 3d hydrogen wavefunctions on each 
atom as a basis set [41]. The exchange hole (x hole) can 
be worked out to be 

^(I*ref j 

thus the exchange hole is actually independent of the po- 
sition of the reference electron. However, the real hole is 
deeper around the reference electron and therefore, de- 
pending on the position of the reference electron, the 
correlation hole (c hole) has to add and remove an equal 
amount of the hole to deepen it around the reference elec- 
tron. Although the xc hole is more localized around the 
reference position, it still shows the two-peak structure 
of the x hole. 

The localization effect becomes more prominent if we 
consider the hydrogen molecule at an elongated bond dis- 
tance of i?H-H — 5-0 Bohr (Fig. 2). Again, the x hole is 
independent of the position of the reference electron and 
is completely delocalized over the molecule. However, 
the real hole is completely localized around the reference 
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FIG. 3. LDA holes for the H2 molecule at -Rh-h = 5.0 Bohr. The reference electron is again at 0.3 Bohr to the left of the right 
nucleus along the bond axis (rrcf = (0, 0, 2.2). Both the LDA x hole and c hole are localized, so do not resemble the exact ones. 
However, their sum, pxc(r|rrof), has much better resemblance to the full xc hole. 



electron, which is again located at 0.3 Bohr to the left 
from the right nucleus. Therefore, the c hole has to com- 
pletely remove the hole from the left side of the molecule 
at put it back at the right side such that it integrates 
still to —1 electron. The c hole can not be regarded as 
a "small" correction to the x hole anymore, since it is 
actually equal in magnitude. The lack of the "small" c 
hole in the Hartree-Fock (HF) description correcting the 
X hole is the main reason for the bad performance of HF 
for the dissociation of molecules. 

It is instructive to consider the LDA holes, since they 
explain why DFT and its predecessor, Xa [42], are so 
successful. The LDA holes are actually the A-integrated 
ones, so a direct comparison is strictly not correct. Luck- 
ily, however, in the dissociation limit the A-integrated and 
the exact hole at A = 1 are identical [8, 43]. The reason is 
that at long bond distances the interaction term A/ri2 for 
A > will favor wavefunction configurations in which the 
electrons are residing on different atoms, i.e. the Heitler- 
London wavefunction. The density corresponding to this 
ground state wavefunction '^x will be exactly equal to the 
wavefunction at full coupling strength. It thus follows 
that '^x = ^' for A > at infinite bond distance. At A = 
the wavefunction simply remains a KS determinant with 
a doubly occupied Ug orbital. However, the A = re- 
gion becomes unimportant [44] in the A-integration for 
the pair-correlation function and hence g = g. Therefore 
at the longer bond distance of -Rh-h = 5.0 Bohr the exact 
xc hole should be close to its A-integrated counterpart. 

In Fig. 3 we show the LDA holes for the hydrogen 
molecule at i?H-H = 5.0 Bohr. From our discussion in 
Sec. II A it is clear what the definition of the xc hole of 
the LDA should be. The following expression [28, 45] is 
consistent with the LDA energy expression (4) 

^(rlr„f) = n(r,ef)(5/i(«(rrcf), |r - r^cfl) - l). (8) 

Gori-Giorgi and Perdew made an accurate model for 
the pair-correlation function of the homogeneous electron 
gas [46] , which we used to calculate p^^^. The x hole can 
be calculated by using the exchange part of the electron 
gas pair-correlation function in this expression and the c 
hole is simply defined as the difference between the other 



two. The most striking feature of the LDA holes in Fig. 3 
is that the x hole is localized instead of delocalized over 
the two atoms as the exact x hole. Although the LDA x 
hole does not resemble the exact x hole at all, the full xc 
hole (the one of interest) is actually modeled quite well. 
Especially, if we consider the exact x hole which is the 
hole employed in HF theory, the LDA xc hole provides 
a large improvement. We also see that the LDA c hole 
only provides a small correction to the x hole: it removes 
the outward oscillations and localizes the LDA hole a bit 
further. Since the correction from the c hole is so small, 
we can understand why the old Xa method already out- 
performed HF so much. In particular, the deepening of 
the hole was empirically taken into account by scaling the 
exchange prefactor, a, from its theoretical value, 2/3, to 
0.7. These observations concerning the LDA holes also 
make it clear that it does not make sense at all to add 
"exact exchange" to LDA correlation. The same holds 
for GGA holes, since they are quite similar to LDA holes, 
only having slightly more wild oscillations and a discon- 
tinuity due to their cut-off to satisfy the sum-rule (3) as 
additional features [47, 48]. 

III. THE SCREENED EXCHANGE MODEL 

Considering these facts about the x hole in dissociat- 
ing H2, it seems to be unwise to add a correlation hole to 
an exact exchange hole. It will be hard to build a model 
for the correlation hole with the proper strongly varying 
features. But we also know that exchange effects often 
dominate and that correlation effects only provide a mod- 
ification of the former. This observation suggests that we 
should not add correlation to exchange but rather modify 
the shape of the x hole by some correlation factor. From 
the holes for the hydrogen molecule (Figs 1 and 2) we 
observe that the main effect of correlation is to localize 
the X hole. This is not special for the H2 molecule, but is 
the main feature of correlation in any system, even in the 
electron gas where the correlation reduces the long range 
behavior of the x hole from to [32]. A further 
example was provided long ago by J.C. Slater when he 
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pointed out that atomic term energies were often accu- 
rately described by term dependent Hartree-Fock theory 
("exact exchange ") by simply reducing Slater's F^. and 
Gk integrals by ~ 25% [49]. 

Following the discussion above it seems rather natural 
to use the following Ansatz for the xc energy E^c of an 
inhomogeneous system 

^>^c = -4/dr|dr'^^f^F(r,r'). (9) 

Here, 7s(r,r') is the spin-integrated non-interacting den- 
sity matrix of the KS orbitals (7) and the "screening" - 
factor F is intended to provide the necessary modification 
of exact exchange which will take care of the effects of 
correlations. The exact expression for F is 



and by multiplying the numerator and denominator by 
n{r) we immediately see that the screening function 
F(r,r') is symmetric in its arguments r and r', cf. (2) 
and (6). 

Our task is thus to find a reasonable model for F and 
we stress again the importance of keeping the symmetry 
of F(r, r') in order to have an ensuing xc potential with 
the correct — 1/r tail outside finite systems. We also 
believe it to be essential to have a model which satisfies 
the sum-rule for the xc hole and in terms of F, our model 
should thus obey 

Jdr' |7,(r,r')pF(r,r') = 2n(r), 

If we, like the founding fathers (KS), first turn to the 
electron gas, we realize that F, and also 7^ for that mat- 
ter, must be spherical functions of only |r — r'|. In the 
spirit of the older WDAs we could thus attempt such an 
Ansatz for our F function. It is, however, important here 
to stress that in the original WDAs it is almost the en- 
tire xc hole which is modeled in this way, whereas in our 
case we just model a modification of the full, in general 
non-spherical x hole. As mentioned previously, we would 
also like to model the _F-function in the case of the dis- 
sociation of II2 . In the dissociation limit of the hydrogen 
molecule the F-function actually takes the form 

FHL(r,r')«/a(r)/a(r') + /b(r)/b(r') 

in terms of two one-point functions fa and located on 
the different hydrogen atoms. This follows from the fact 
that in the limit of large separation between the nuclei, 
the Heitler-London wavefunction becomes exact — but we 
will not present the details here. But it means that in 
this limit, the i^-function is very far from spherical. 

When the two electrons are on different nuclei, the F- 
factor should vanish, because we are then dealing with 
two non-interacting subsystems and there is neither ex- 
change nor correlation. When both electrons are on 



the same atom, the F-funcion should actually be 2 to 
make the xc hole equal to the negative of the local den- 
sity. In this way the xc hole will precisely remove the 
full self-interaction on each atom — not just half the self- 
interaction as Hartree-Fock does — and we recover the 
correct limit of two separate and non-interacting atoms. 

As suggested by the above observations, the following 
Ansatz for the screening function F{r,r') might stand a 
chance of carrying us from the limit of a homogeneous 
system to that of the complete breaking up of the H2 
bond 

F^^{r, r') = A(r)A(r') h{\r - r'|, n(r, r')) . (10) 

Here, the spherical function h has an effective screening 
radius fg given by Airf^ = 3/n. This form was also in- 
spired by the success of a wavefunction for the helium 
atom by Hylleraas [50] having precisely this form. 

We are then left with the choice of satisfying the hole 
sum- rule either by adjusting the one-point function A(r), 
the hole-depth function, or by varying the effective den- 
sity n{r,r'). We stress that the latter must be a sym- 
metric two-point function in order not to jeopardize the 
asymptotics of the potential. In terms of the Ansatz (10) 
the sum-rule reads 

A{v) J dr' |7,(r, r')pA(r') h{\v - r'|, n(r, r')) = 2n(r). 

(11) 

The effective density is expected to vary in a similar way 
as the density itself (from very small to very large values) 
and it is a two-point function. The hole-depth function 
A on the other hand is a one-point function of limited 
variation — typically between zero and two depending on 
the normalization of the function h. Consequently, it ap- 
pears to be numerically more stable to use the A-function 
for the purpose of satisfying the sum-rule. Indeed, we 
have encountered no difficulties in solving (11) in our ap- 
plications. A further argument in favor of this choice is 
a lack of guidance from the electron gas when determin- 
ing the A-factor, should we have chosen to satisfy the 
sum-rule by varying the effective density n. Most WDA 
models have used the latter procedure, but again, they 
have not considered an A-factor. 

In our case we are thus free to choose the effective 
density. Typical choices are n{r,r') = ^(^(r) -|-n(r')) or 
n{r,r') = n((r-|-r')/2). in previous attempts to construct 
a symmetric version of an WDA-like model we have found 
that the first of these suggestions gave rise to numerical 
difficulties, whereas Gunnarsson et al. [28] encountered 
difficulties with the second choice. A choice which seems 
to work in our previous attempts is 

n{r,r') = ^nir)n{r'), (12) 

which is the definition of the effective density n which 
we will use here. We stress, however, that there is no 
compelling reason for this choice. As a matter of fact. 
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this is one part of our new model which we might have 
to revise in future attempts to improve the accuracy of 
the modeL 

It is not clear what properties and shape the screening 
function h should have. However, since its main task is 
to screen the x hole, we will use the very simple form 
inspired by the screened Coulomb (Yukawa) interaction 

/iHEG(r,n)=e-^(«)''". (13) 

The function D{n) is fitted to the homogeneous electron 
gas such that it yields the exact xc energies for all densi- 
ties. More details on the fit can be found in appendix B. 
This form for the screening function is definitely too sim- 
plistic. More knowledge is required to build more accu- 
rate SX models and will be the subject of future research. 



IV. RESULTS 

One of the most severe tests for our new SX functional 
will be if it performs well for dissociating molecules. To 
keep matters simple, we have limited ourselves to the 
hydrogen molecule. The most natural grid to perform 
calculations on a diatomic molecule is an elliptic grid 
where the foci are at the positions of the nuclei and then 
to rotate the elliptic grid around the axis connecting the 
two foci, the z-axis, i.e. a prolate spheroidal coordinate 
system. More details on the used grid can be found in 
appendix C. As a further simplification we used the den- 
sity from a full CI calculation in the same basis as used 
before (Is, 2s, 3s, 2p and 3d hydrogen wavefunctions). 
This CI expansion is not very good for obtaining an ac- 
curate total energy. However, we believe that the density 
will be accurate enough for the SX model. 

The first step in evaluating the SX model is to solve 
for the hole-depth function, A{r), by solving the integral 
equation (11) on the grid. Once the hole-depth function 
is obtained, we performed the double integral (9) with 
F^^ (10) to obtain the xc energy according to our sim- 
ple SX model. To calculate the total energy, we note 
that the non-interacting kinetic energy can be directly 
obtained from the density for the two-electron system, 
T5 = I / (V-\/n) , and the Hartree and nuclear contribu- 
tions are already known from the full CI calculation. In 
Fig. 4 we compare the total energies from the SX model 
with the ones from a full CI calculation as a function 
of the distance, i?H-H, between the hydrogen atoms. As 
mentioned before the Slater basis is too poor to obtain a 
good total energy, so we performed a full CI calculation 
with the DALTON package [51] in an aug-cc-pVQZ ba- 
sis [52] as a reference and additionally the corresponding 
HF result is shown as well. Unfortunately our SX model 
is not performing much better than the HF functional. 
It follows the HF curve rather closely and only around 
^^H-H ~ 6 Bohr the curve starts to bend downwards to 
the correct total energy. The most important feature of 
the SX model is that it directly models the xc hole and 




1 23456789 10 
Rh-h [a-u.] 

FIG. 4. Comparison of the total energy from the simple SX 
model with the ones from the full CI calculation for varying 
bond-length. 



therefore we can look at it what is going wrong. In Fig. 5 
in the left panel, we show the hole at i?H-H = 1.4 Bohr. 
If we compare with the exact holes in Fig. 1, we see im- 
mediately that our current SX model is not localizing 
the X hole sufficiently. Actually, the SX hole is almost 
identical to the x hole px ■ In the left panel of Fig. 6 we 
show the hole the elongated distance i?H-H = 5.0 Bohr. 
Comparing with the exact holes in Fig. 2 we see that 
the SX model actually does localize the x hole, but not 
sufhciently. The lack of localization of the hole explains 
exactly why the total energy in the SX model is consis- 
tently too high (Fig. 4). Due to the l/|r — r'| term in the 
expression for the xc energy (1), a too delocalized hole 
does not stabilize the energy enough, which results in a 
too high energy. 

In retrospect we should not be surprised that 
parametrizing the screening function by the homoge- 
neous electron gas did not work out. The main task 
of the screening function is to localize the hole which is 
most important in inhomogeneous systems like the hy- 
drogen molecule. Therefore, its actual form and localiza- 
tion strength should be modeled by these systems and 
not the homogeneous electron gas. A detailed study and 
parametrization are beyond the aims of this article, but 
to reinforce our arguments for this statement, we also 
like to show some results for the following two heuristic 
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h„„ h, hj LDA exact (A=l) 




FIG. 5. The different model holes of the H2 molecule at -Rh-h = 1-4 Bohr and the reference electron at 0.3 Bohr to the left 
of the right nucleus along the bond axis (rref = (0,0,0.4) Bohr). The positions of the nuclei are indicated by the blue lines 
and the position of the reference electron is indicated by the red line. The right most panel shows the exact (not integrated) 
xc hole for comparison. From left to right, the panels show the xc holes from the SX model with the /iheg screening function 
(13), the SX model with the hi screening function (14a), the SX model with the h2 screening function (14b) and the LDA hole 
evaluated as in (8) with the pair distribution from Ref. [46]. 



h„„ h, hj LDA exact (A=l) 




FIG. 6. Similar to the previous plots, but now for -Rh-h = 5.0 Bohr. The reference electron is still at 0.3 Bohr to the left of 
the right nucleus along the bond axis (rref = (0,0,2.2) Bohr now). 



screening functions 

hi{ri2,n) ^ exp(^-ci(J^^^, (14a) 

ft2(ri2,n) =exp(^-C2(y^)'y (14b) 

The first one, hi, keeps the Slater like form, but replaces 
the parametrization by the electron gas by a simple di- 
vision by fs to make the total dimensionless and a con- 
stant ci that we can choose to our liking. We found that 
ci = 2.0 gave a nice dissociation behavior for the energy. 
The second one is mainly included to emphasize that the 
required shape of the screening function is also unclear 
at the moment. Choosing the constant in the same man- 
ner as before, we found C2 = 0.5 to be sufficient for our 
purposes. Note that both screening functions are not fit- 
ted to the electron gas anymore and are therefore not 
expected to give the correct xc energy density, Exci^s), 
for the gas. 

The results for the energy from these screening func- 
tions are shown also in Fig. 4. The situation is now 
rather different. The total energy is consistently underes- 
timated. However, the shape of the curve is definitely an 
improvement. The total energy at elongated distances 
_Rh-h > 5 Bohr remains rather constant as we choose 
the constants Ci to do so. From the figure it is not im- 
mediately clear if the shape is also an improvement over 
the simple LDA functional. However, shifting the curves 
such that their minima coincide with the full CI mini- 
mum (Fig. 7), we see that the energy from the SX mod- 
els follow the full CI energy much closer. In particular 
the Gaussian, /12 seems to reduce the overbinding most. 
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FIG. 7. Comparison of the total energy of the LDA func- 
tional [46] and the SX model with the Gaussian, /12, as its 
screening function with the full CI results. The total energy 
of the LDA and SX model are shifted such that the minima 
coincide with the full CI value. 



Even at equilibrium distance the position and shape of 
the minimum seems to be somewhat improved. 

Again, since we have built a direct model for the hole 
we can explain both features. Considering the holes from 



9 



the heuristic hi and /12 screening functions in Figs 5 
and 6, we see that both screening function are too power- 
ful: they screen the x hole too much. This explains why 
the total energy is consistently lower than the full CI re- 
sult: the hole becomes too compact. Since the hole is 
even more compact for hi screening function than the h2 
screening function, the energy of the hi screening func- 
tion is even lower than the energy of the /i2 screening 
function. However, if we consider the shapes of the holes, 
then they are definitely an improvement over the previ- 
ous screening function. Especially at the elongated bond 
distance Rh-h = 5.0 Bohr, we see that both heuristic 
screening functions nicely localize the hole completely at 
the correct side, which explains their improved energy in 
the dissociation limit over the erroneous l/i?H-H of HF. 

The LDA hole is also shown in Figs 5 and 6. Compared 
to the HF, cf. px in Fig. 1 and 2, the LDA hole is a big 
improvement since it localizes around the reference elec- 
tron. This improvement is also visible in the total energy 
where the LDA does not exhibit the erroneous l/i?H-H 
behavior as HF does (Fig. 4). However, the shape of the 
LDA hole is ridiculous. It is spherical by construction 
and it does not have the peaked features at the nuclei. 
Instead, the LDA hole attains its minimum at the posi- 
tion of the reference electron. The SX holes, especially 
with the /12 screening function, are an improvement over 
both the HF hole and the LDA hole. It correctly localizes 
the hole around the reference electron and still retains the 
distinct peaked features of the hole. Hence, the binding 
curve is much improved over HF and LDA (see Fig. 7). 



V. CONCLUSION 

The aim of the present work is to construct a func- 
tional for the exchange-correlation (xc) energy of DFT 
which applies to such diverse systems as the electron gas 
and the dissociation of the hydrogen molecule (H2). To 
this end we try to find a model for the xc hole of these and 
intermediate systems in real space. It has been known for 
long that the exchange (x) hole captures many important 
features of the full one. For instance, in atoms term en- 
ergies are often well described by reducing the exchange 
effects by 25% and in the electron gas correlations have 
the effect of reducing the range of the pure x hole from 
a decay to a decay, r being the distance from 
the the center of the hole. Thus unlike previous models 
that have sought to model the xc hole in real space, our 
present model aims at modifying the full x hole. Our 
investigations have shown that in the case of the disso- 
ciation of H2 the xc hole is qualitatively different from 
the X hole. Therefore, it is an unwise strategy to add a 
correlation (c) hole to a full x hole. The former would 
have to replace the latter with a full xc hole with appro- 
priate features. We show here that this can be achieved 
in a more natural manner by multiplying the x hole by an 
appropriate correlation factor. Judging from the electron 
gas, the correlation factor might be chosen as a function 



of the distance |r — rj-of | between an electron and a refer- 
ence electron. (We here again remind the reader that the 
xc hole describes the depletion of negative charge around 
an electron known to be at the reference position ricf.) 
Unfortunately, our investigations have shown that such 
a simple correlation factor will have difficulties in mov- 
ing half an x hole from one atom to the other — which is 
the appropriate effect of correlation in H2 at large nu- 
clear separation. Instead, we have chosen to include in 
our correlation factor an additional factor A{r) depend- 
ing on only one coordinate — a modulation of the depth 
of the xc hole. Such a factor is by symmetry just a con- 
stant in the electron gas and irrelevant to the shape of 
the hole in that case. Consequently, we have no guid- 
ance from the gas in choosing the A-factor. Instead, we 
have decided to determine this factor at every point in 
space from the sum-rule for the xc hole. This sum-rule 
is of course of utmost importance for obtaining an xc 
potential with the correct asymptotics outside finite sys- 
tems (— 1/r) from our model. This important property 
also requires the full symmetry in r and r' in the density 
multiplied xc hole, n(r)pxc(r'|r), something that we have 
emphasized throughout the paper. 

The decision to use the A-factor for satisfying the sum- 
rule leaves us with great freedom in choosing a screening 
factor depending only on |r — r'|. Thinking about the 
electron gas it is natural to let this screening function 
have a screening length fg determined by an effective 
density n according to Annf^ = 3. In the present work 
we have made the somewhat arbitrary choice n(r, r') = 
V^n(r)n(r'). We are, however, aware of that we might be 
forced to abandon this simple choice in future refinements 
of our model. 

For the actual shape of the screening function we have 
simply made a couple of reasonable choices for illustra- 
tional purposes. They are rapidly decaying, analytic 
functions with a few parameters with density depen- 
dence. One of the screening functions had its parame- 
ters specifically chosen such to reproduce the "exact" xc 
energies of the electron gas in the homogeneous limit, 
whereas two other screening functions had more ad-hoc 
parameters to illustrate the effect of selecting different 
forms of the screening function. 

We could, however, also have attempted to find a 
screening function with a shape that would have allowed 
us to obtain accurate results for one- and two-electron 
systems. We would then have had an approximation 
which is able to properly dissociate H2, which would be 
exact for the electron gas and quite accurate for the one- 
and two-electron cases. Such a functional is likely to give 
good total energies in a large number of systems. 

In the present work we have, however, refrained from 
going down this road. Instead, our aim here was to 
present the basic ideas and to leave further refinements to 
future investigations. In order to just illustrate our new 
approach, we thus chose to present results obtained from 
two simple, but rather arbitrary screening functions given 
by the equations (14a) and (14b). We have seen that the 
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shorter ranged choice (/i2) gives better results, but they 
are still not very good. It is however, seen from both 
Fig. 4 and Fig. 7 that the errors are mainly associated 
with an inadequate description of a simple one-electron 
system. (At a bond distance of 10 Bohr we basically have 
two separate hydrogen atoms.) 

The most important result of the present work is that 
we managed to design a model which is able to describe 
the proper behavior of the xc hole of a hydrogen molecule 
as it dissociates. The details are not overly accurate, but 
we nurture real hope that they may fall in place by a 
better choice of the screening function. For now, however, 
we have left the search for such a function to future work. 



Appendix A: The xc potential 



differentiation. Working out the first part of the poten- 
tial gives 

SO we find the correct Coulombic 1/r behavior. We 
only have to check the charge. Using 75(ri,r2) = 
n(ri)n(r2) in the sum-rule for the hole-depth func- 
tion (11), we find 



A(r) 



J dr' 72{r')A{r')h{\r ~ r'|, n) = 1, (A3) 



so indeed, the Coulombic part of the potential also has 
the correct charge of —1. 



The xc potential is obtained by straightforward differ- 
entiation of with respect to the density 
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/ ^l7.(r.r2)P 
\ dn{r) 

+ 2|7,(ri,r2)P^^A(r2)%i2,n) 
dn(r) 

+ |7.(ri,r2)r^(rOA(r2)^^^|;|^). (Al) 

The functional derivative of the hole-depth function can 
be obtained by differentiating its hole sum-rule (11) with 
respect to the density. Its derivative can be worked out 
as 



M(ri) , A2(ri) /•, , , ,,2, , _,M(r2) 
dn[r) 2n(ri) J dn(r) 



Air) 



n{r) 
I7s(ri,r2)p 
Sn{r) 



S{r - ri) - 



2n(ri) 



i| y"dr2 A{V2) X 



ft.(ri2,n) + |7s(ri,r2)| 



2(5/i(ri2,n) 



Sn{r) 



Unfortunately this equation does not give a closed expres- 
sion for the functional derivative of ^(r), but just like its 
original counterpart (11) has to be solved iteratively to 
self-consistency. 

Since the hole correctly integrates to —1 electron (11), 
the potential has the correct asymptotic behavior — l/r. 
However, due to all the functional derivatives this is not 
directly visible in the expression for the xc potential (Al). 
The main complication arises from the functional deriva- 
tive of the KS density matrix whose evaluation requires 
the application of the chain-rule multiple times. Fortu- 
nately, in the case of a two-electron system, the Kohn- 
Sham density matrix can be expressed directly in the 
density as 7s(ri,r2) = y^n(ri)n(r2) allowing for direct 



Appendix B: The electron gas screening function 

One of the requirements of the new functional is that 
it should work form a broad class of systems. Not only 
for molecules, but also for extended systems and in par- 
ticular the homogeneous electron gas. Since the homoge- 
neous electron gas is well studied and much of its prop- 
erties are known, it provides the ideal system to fit the 
screening function (13) such that the SX functional will 
be exact for the homogeneous electron gas, in particular, 
it should give the exact xc energy for the electron gas. 
Before we can start fitting the screening function, we 
have to solve for the hole-depth, A. For the KS density 
matrix of the homogeneous electron gas one can straight- 
forwardly calculate that 



kp sin(fcFri2) - fcFf'i2 cos(fcFri2) 



{kpru 



where the Fermi wavevector kp 



Sir^n. Using this ex- 



pression for the KS density matrix in the sum-rule for the 
hole-depth function (11) and the Ansatz for the screening 
function (13), we find 



1 = -A'F^iD), 

TT 



(Bl) 



where we defined D = D/kp and the integrals 



F„(/3)= / dy(sin(j;) 



y cos 



(y))\~ 



More details about these functions and explicit expres- 
sions are given in Appendix D. Now we have an explicit 
relation how to calculate the hole-depth function A for 
the homogeneous electron gas (Bl), we will fix D by re- 
quiring our functional to produce the exact xc energy for 
the electron gas. In terms of our model, the xc energy 
density becomes 



-—A^F^iD). 



(B2) 
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Since the equations (Bl) and (B2) are linear in A^, the 
hole-depth function can simply be eliminated by dividing 
the equations, which gives 



F5{D) 3 £,c(r.) 



(B3) 



FiiD) 27r e^{rs) ' 

where we used the Seitz radius = 3/(47rn) and that 
Ex = —Skp/i^Tr) for the homogeneous electron gas. To 
find D{rs), we need an explicit expression for £xc(^s) for 
which we used an accurate fit to the random-phase ap- 
proximation (RPA) and Green's function Monte Carlo 
data by Perdew and Wang [53]. Unfortunately, the ex- 
pression for D{rg) cannot be inverted analytically. How- 
ever, one can show that Fc,[ji) / Fi{ji) is a monotonic in- 
creasing function over /3 > 0, so at least the solution 
to (B3) is unique. Further, in the low density limit we 
have 



,. 3 ey,c{rs) 
hm — - 

r.^oo 27r Exit's) 




0.92925, 



where the parameters ai = 0.21370 and (3^ = 0.49294 
are from the low-density fit of Perdew and Wang [53]. 
Inverting (B3) numerically, we find that in the low den- 
sity limit 



D{rs) 



where D°° — 2.27591. For the asymptotic behavior in 
the high density limit, — !■ 0, we find for the ratio of 
the xc and x energies 



3 £xc(^s) 



3 

2tt 



^ i:i 2l 



-Co rs \ii{rs) + 



27r e^{rs 

where cq = Cq^^ = (l — ln(2)) /tt^ is a constant from the 
high-density fit by Perdew and Wang to the RPA [53]. 
Similarly for small arguments of -F5/-F4 we have 



3 

2tt 
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27r2 



Comparing these two high density limits, we find that for 
high density to lowest order 



^9(l-ln(2)). 



To obtain a workable expression for D{ra), we simply 
solved (B3) numerically for several Tg and made a least 
square fit using the following the following Pade approx- 
imant 



ao + airs + b^D^^r^ 



for which we found the coefficients oq = 0.149056, 
fli = 0.180374, 61 = 1.16435, 62 = 0.128538 and 63 = 
0.000703698. Note that we did not enforce the proper 
limit for Ts — >■ 0, since the low density region is more 
important and relaxing this constraint significantly in- 
creased the accuracy for Vg > 0.1 Bohr, which are more 
relevant densities in non-relativistic molecules and solids. 



Appendix C: The prolate spheroidal grid 

In this appendix we very briefiy introduce the prolate 
spheroidal grid and give details on the exact parame- 
ters used for the grid in the calculations. The prolate 
spheroidal grid is created from an elliptic grid by rotat- 
ing it around the axis connecting the two foci, the z-axis. 
The elliptic coordinates are defined as 



ri + 7-2 
2p 



and 



77 



ri - r2 
2p ■ 



where is the distance to nucleus i, which are placed 
at ±p from the origin on the z-axis. One readily sees 
that the ranges of the elliptic coordinates are 1 < C £^nd 
~1 < < 1- The curves ^ — constant describe ellipses 
and the 77 — constant curves describe hyperbolae. The 
intersection points of the ellipses and hyperbolae for dif- 
ferent constants will be used as grid points. Also taking 
the revolution around the z-axis into account one can de- 
rive expressions for the cartesian coordinates in therms 
of the prolate spheroidal ones 



V = p^/{e-m-v^)^H<^>). 



p^-q. 



The main advantage in using ellipses and hyperbolae is 
that they are orthogonal to each other. Therefore, the 
metric g is simply diagonal, ~ g^^j, = g^p^ = 0, with 
diagonal elements 



5« 

9nn 



dr 

dr 

drj 
dr 

84) 



dr\ 
dr \ 



.2^ 



1 ' 

^2 



1 



The scaling factors are defined as Ai = y/gu, from which 
we immediately find the volume element 

and the gradient 



V = 



]_d_ 

\x^ dcpj 



( 



e - 772 d( 

/ 1 ~ 77^ d 

1 



VV(e^-i)(i-772)a0/ 



We could also write down an explicit expression for the 
Laplacian, but we do not need it. The disadvantage of 
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discretizing the Laplacian directly is that it is not sym- 
metric anymore. Instead we use that for functions van- 
ishing sufficiently fast at the boundary [54] 

- y"dr/(r)V^g(r) = ^dr V/(r) • V.g(r), 

so formally we can write 

which is inherently symmetric and its discretization can 
directly be constructed from the discretized gradient. 

The ^-grid points are generated by starting from an 
equidistant grid, u € [0, 1), and using a simplified version 
of the transformation used by Becke [54] to generate grid 
points and weights for the ^-grid 



an) = ^ 



1 



(1 - u2)P5 ■ 

The parameter can be used to affect the distribution of 
the grid points between the inner and outer region. For a 
given maximum value of ^, ^max, is easily determined 

as 



Pi 



ln(Aw(2 - Aw)) ' 



where Au is the distance between the points in the u- 
grid. An important aspect is the which has the benefit 
that the cusps at the nuclei are transformed into smooth 
Gaussians in it-grid. The disadvantage is that the weight 
becomes zero at the line between the nuclei {w{£, = 1) = 
^'{u = 0) = 0), so the solution cannot be calculated 
directly at these points and has to be extrapolated. 

For the ?7-grid we started from an equidistant grid, v e 
[0,7r], and use the following transformation to generate 
the ?y-grid points 



ri{v) 



cos 



V + p,f sin(2t;)). 



where the parameter can be used to modify the point 
distribution between the intermediate region and the nu- 
clei. Becke found Prj = —0.25 to be optimal [54]. Also 
this transformation has the property that the cusps at 
the nuclei are transformed into Gaussians in the u-grid. 



As one might already expect, the disadvantage of this 
feature is that the weights at the boundaries vanish 
{'w{ri = —1) = w{ri = 1) =0), so the results will have to 
be extrapolated also to these points. The rj — ±1 points 
correspond along the bond axis outside the molecule. 

We found that 80, 81 and 40 points in the ^, rj and 
^-directions respectively gave numerically sufficient con- 
verged results. For the self-consistent solution of the 
hole-depth equation from its sum-rule (11) it was impor- 
tant that the density did not become too small. Hence 
the grid should not have points far in the asymptotic re- 
gion. Using ^max = 1 + 10/p was sufficient to prevent 
points with numerically vanishing density, though still 
including a sufficient part of the relevant space. 

In the equidistant grids it is straightforward to define 
numerical derivatives. For example one can use cubic B- 
splines [54] or simple central finite differences [55] which 
we used in our calculations. In practical calculations, 4th 
order derivatives already gave sufficient accuracy. 



Appendix D: The functions Fn 

In this appendix we show how the integrals defining 
the functions F„ can be evaluated. The functions F„(/3) 
satisfy 



dj^n+i(/3) 
d/3 



(Dl) 



Using the boundary condition _F'„(+oo) ~ 0, we find 



(D2) 



which can be used to obtain successively higher order 
functions. The integral most suitable for direct evalu- 
ation is Fq, albeit the evaluation is still rather tedious. 
The final result is 



5fl2 - 



4 



(D3) 



By applying successively the integral formula (D2) we 
can obtain the higher order integrals required for our 
model. Their evaluation is straightforward, but takes 
the necessary amount of time. The final results are 
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ln(l + 4//32) 1 4 3/2 



F3W) 

F5W) 
FeiP) 



2/32 (/32 + 4)2 /32+4' 



i^^2(/3) = ^ arctan(/3/2) ' ^ ^ 



4 '"^^ ' 2/32 +4' 

(/32 + 2) ln(l + 4//32) _ 1 
8 2' 
arctan(2/^) (/^^ + 6/3) ln(l + 4//32) ^ 

3 24 ^ 6 ' 

1 /32(/32 + 12)ln(l + 4/^2) ^2 /3arctan(2/^) 

4 96 24 3 ' 

(5/32 +4)arctan(2//3) ;93(-^2 _^ 20) ln(l + 4//32) 11/3 
30 480 ^ 120 ^ "60"' 



The integrals for n > 6 do not converge anymore, so full CI program. K.J.H.G. and R.v.L. acknowledge the 
Fq{/3) is actually the last function in this series. Academy of Finland for research funding under Grant 



No. 127739. 
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